Fractal iso-contours of passive scalar in smooth random flows 
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We consider a passive scalar field under the action of pumping, diffusion and advection by a smooth 
flow with a Lagrangian chaos. We present theoretical arguments showing that scalar statistics is 
not conformal invariant and formulate new effective semi-analytic algorithm to model the scalar 
turbulence. We then carry massive numerics of passive scalar turbulence with the focus on the 
statistics of nodal lines. The distribution of contours over sizes and perimeters is shown to depend 
neither on the flow realization nor on the resolution (diffusion) scale ra for scales exceeding r^. 
The scalar isolines are found fractal/smooth at the scales larger/smaller than the pumping scale 
L. We characterize the statistics of bending of a long isoline by the driving function of the Lowner 
map, show that it behaves like diffusion with the diffusivity independent of resolution yet, most 
surprisingly, dependent on the velocity realization and the time of scalar evolution. 

PACS numbers: 47.27.-i, 47.51. +a, 47.27. wj 



A fundamental problem in mixing is to describe the 
iso-contours of the quantity which is mixed (called pas- 
sive scalar in what follows). This is important for nu- 
merous practical applications. Renewed interest in this 
problem is related to the recent mathematical advances 
in describing random curves, particularly identifying an 
important class of curves SLE K that can be mapped into 
a one-dimensional Brownian walk with the diffusivity k. 
SLE stands for Schramm-Lowner Evolution and presents 
a fascinating subject where physics meets geometry [J- 
3]. If lines belong to an SLE class, as was found, in 
particular, for isolines of vorticity and other quantities 
in turbulent inverse cascades [J-|6|, conformal invariance 
provides exact formulae for the statistical description of 
the kind unimaginable before in turbulence studies. 

It may seem that the simplest setting for mixing is 
when the velocity is spatially smooth, which one encoun- 
ters in many natural and industrial flows. An important 
problem of that type is large-scale mixing in the Earth 
atmosphere where the turbulence spectrum is approx- 
imately as steep as k~ 3 0); the velocity gradients arc 
then well-defined and the flow can be locally considered 
smooth. The same type of flows one encounters in the 
phase space of dynamical systems. Mixing in smooth 
flows is provided by exponential separation of trajecto- 
ries and Lagrangian chaos. We consider an incompress- 
ible fluid flows which correspond to Hamiltonian flows in 
phase space. Generally, such flows have a positive Lya- 
punov exponent A so they stretch any element into a long 
narrow strip. Additionally, our passive scalar is subject 
to random sources/sinks distributed in space with a cor- 
relation scale L. Resulting "passive scalar turbulence" 
in this (so-called Batchelor [§]) regime has been exten- 
sively studied in terms of the one-point probability and 
multi-point correlation functions [9(, one can even write 
down a closed expression for the probability of any given 
scalar field realization. Yet to the best of our knowledge, 



nothing is known about the properties of an infinite-point 
object, isoline. This may seem surprising since much is 
known about isolines in a non-smooth (fully turbulent) 
velocity field [1, Il0l - fl2| . Apparently, the description of 
contours in the Batchelor regime is more involved. A con- 
ceptual difficulty is related to the lack of scale invariance. 
Indeed, every long contour is simultaneously stretched to 
the scales far exceeding L and contracted to the scales 
much less than L in transverse directions. Technically, 
experimental and numerical studies of passive scalar tur- 
bulence in the Batchelor regime are notoriously difficult 
because of a slow logarithmic convergence of the correla- 
tion functions with the resolution 1 1 31 ] . 

In this work we suggest a new, very effective, method 
of numerical simulation of the passive scalar turbulence 
in the Batchelor regime. The method is based on the an- 
alytic representation of the Lagrangian path integral. We 
carry out extensive numerical simulations of the passive 
scalar turbulence in two dimensions and find out that the 
isolines are non-fractal one-dimensional lines at the scales 
less than L. We then show that the isolines are fractal 
at the scales larger than L, and describe the statistics 
of the contour sizes and perimeters. Finally we explore 
relations of these isolines to SLE. 

Passive scalar 9 is carried by the velocity v, forced by 
ip and diffuses with molecular diffusivity k^: 



d t + (t> • V)0 = K d V 2 d + ip . 



(1) 



We take the pumping ip(r,t) to be white Gaussian with 
zero mean and variance tp(r, t)ip(0, 0) = x(r)<J(f), where 
X(r) decays faster than any power at r > L. Our new 
computational method of generating the scalar field ex- 
ploits linearity of the advection-diffusion equation (flj). 
To get a snapshot of 0(r,t), we sum over all of the 
blobs of scalar that hit the surface at random positions 
and random times in the past (v times per unit time), 
sec Figure [T] Each blob at its initial time to has an 
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FIG. 1. A snapshot of a part of the passive scalar field of size 
(75 L) 2 , where L is the forcing scale. The color bar shows the 
magnitude of 9. 



isotropic shape centered around a random position r c 
and of random amplitude Qq: 0(io,?"o) = ©oexp[— (to — 
r c ) 2 /(2L 2 )]. The shape of such a blob at time t is found 
by characteristics: 



r ) = 



e i 2 



det I(t,t ) 



=e 2 



i(r(t)-r )/- 1 (t,i )(r(i)-r (! ) 



, (2) 



where r(t) = W(t, to)r(to) defines the evolution operator 
W(t, to): obtained for a white in time and spatially linear 
velocity profile with 



dVjjt) dv k {t>) 
dxj dxi 



DS(t' - t)[35ijSji - SijSki - SuSkj] 



The moment of inertia is I(t,t Q ) = W(t, t )W(t, t ) T 

+n d /* dt'w(t,t )w(t' ^or^w^toW- 1 ^' ,t )} T '. 

Seven values specify a blob at time t: the symmetric 
matrix I(t,to), 8o, r c and <o- The blob database is 
made large enough to ensure Gaussianity of 9. To obtain 
the isolines of 9 we used MATLAB's contourc function, 
an example of the contour can be seen in see Figure [2] 
We measure the box-counting generalized fractal dimen- 
sions: D q = lim e _K) logEf {e) P q i{e))[{q ~ !) log(e)] -1 - 
Here Pi{s) is the probability of finding a point in the 
i— th square of area e 2 , and N(e) is the number of e 2 
squares needed to cover the contour. For long contours, 
we study the statistics of bendin g by making the Lbwncr 
map of the curve into a line pi Il4|. As one goes along 
the curve, the image moves along the line; if this motion 
(called driving function) is Brownian then the curve 
belongs to SLE. We made several velocity realizations 
with different resolutions r d /L e {0.06,0.03,0.015}, 
pumping frequencies vj\ e {0.05,0.02,0.01,0.005} and 
scalar evolution times TX € {20, 50, 100}. 

Let us now confront theoretical expectations with the 
results of numerics. One expects that pumping alone 




FIG. 2. A long scalar iso-contour in a window 24001/ x 360L, 
where L is the forcing scale. 



would produce a Gaussian field 9 whose zero isolines are 
smooth at the scales below L, while at larger scales they 
are equivalent to critical percolation called SLEg (hav- 
ing the dimensionality D = 7/4). Velocity field by itself 
does not change the statistics of 9 as it just rearranges it; 
the flow stretches isolines uniformly at the direction of 
the eigen- vector of the positive Lyapunov exponent and 
contracts them transversal to it. Non-trivial statistics 
of 9 and its isolines arises from an interplay of velocity, 
pumping and finite diffusivity or finite resolution, which 
leads to the dissipation of 9 and reconnection of isolines 
that came closer than the resolution scale r d - We as- 
sume r,j < i. At the scales between L and r ( i, there 
is a cascade of passive scalar whose correlation functions 
of the scalar are logarithmic: (9 n (0)9 n (r)) - \n n (L/r) 
0. Lower orders, n < ln(L/r), correspond to Gaussian 
probability density function (PDF), while the PDF tails 
are exponential Q. The scalar field itself is thus non- 
smooth at r < L, what about its isolines? If the scalar 
was a Gaussian (free) field with logarithmic correlation 
functions, it would have the isoline with the fractal di- 
mensionality 3/2. Let us show that the Gaussian PDF, 



V{9} oc exp 



-| / dr|V0(r)| 2 /> 



(3) 



does not satisfy the respective Fokker-Planck equation: 

6 2 



dt 2jJ Vl V2 69( ri )S9(r 2 ) 



[K a/3 V a 9( ri )Vp9(r 2 ) 



-X(ri - r 2 )]V + n d J dr 1 ^^A9(r 1 )V. (4) 

Indeed by substituting V in @ we see that already the 
first term gives a non- vanishing contribution (highest or- 
der in 9) 

JJdn dr 2 DV{2[V9( ri ) ■ V0(r 2 )] 2 

-|V0(n)| 2 |V0(r 2 )| 2 }. (5) 

Indeed, we know that correlation functions of 9 include 
cumulants (l5l |. One may argue that the cumulants con- 
tain less logarithmic factors than reducible terms and are 
small @. However, the properties of isolines must de- 
pend on those cumulants, since they contribute the cor- 
relation functions of the gradients in the main order. It 
is straightforward to establish that the correlation func- 
tions of the passive scalar are not conformal invariant, 
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i.e., for instance, the four-point function does not have 
the form 

Fa = J [ , (r^r^ri^zri^A) . (6) 

Indeed, we know from [l5| that F± = F(r 12 , r^) + 
F{ r i37 r 2i) + F{ r i4, r 23) which is compatible with © 
only for Gaussian statistics, which is not the case, as we 
have shown. Therefore, passive scalar is not in any way 
close to a free field. 

Note in passing that one can also show that the 
scalar statistics is not conformal invariant in a com- 
pressible Id Batchelor-Kraichnan model, where veloc- 
ity is Gaussian with the zero mean and the variance 



(v(x,t)v(0,0)) = 5(t) 



K 



Correlation function 



of the scalar satisfy (d t + J2 x ij^i^j)(^^2) = xi x i2)- 
At 1 > L the pair correlation function of the scalar 
has parts linearly growing with time while the struc- 
ture functions are finite: {{8\ — O2) 2 ) ~ x(0) ln(xi2/£) 
fl6j ] . The two-point correlation function of the gradients 
is (u>iU>2) ~ x(0) x i2- The four-point correlation function 
of the gradients can be calculated exactly for x( x ) K 
cxp(— x/L) using the method of [jq : (CJ1W2W3W4) = 



x 2 (o)£[; 



mn^kl 



+ 2(x mn + X k l) 2 /(x mn Xki)} for TO jt 



n 7^ k 7^ I. Both parts (Gaussian and cumulant) are zero 
modes of the operator C = ^didjX 2 ^. The cumulant 
is not conformal invariant: under the transformation of 
coordinates x' = (ax + b)/cx + d) with ad — be = 1 it 
changes its form rather than just acquire the factors (Ja- 
cobians) (dx' /dx) A = l/(cx + d) 2A . The dimension of uo 
is A = 1 as can be seen from the invariance of the pair 
correlation function: (xi—X2)~ 2 {cx\+d) 2A (cx2 + d) 2A = 
(x[ - x' 2 y 2 = (xi - x 2 )~ 2 {cxi + d) 2 (cx 2 + d) 2 . We thus 
conclude that the scalar field is not conformal invariant. 

If one tries to find an analogy with a non-Gaussian 
field having logarithmic correlation functions, such is the 
height function built on independently oriented loops 
from the 0[n] model [l7j |. deviations of |« — 4| and 
\D — 3/2 1 are proportional to cumulants in this case. 
Another much exploited similarity is between the pas- 
sive scalar and the vorticity cascade of two-dimensional 
turbulence; vorticity isolines was shown to have multi- 
fractal isolines with dimensionalities changing from 3/2 
to 1 ||. Despite all these suggestive similarities, let us 
show that, contrary to these expectations, scalar isolines 
are smooth below L. Figure [3] shows the box-counting 
dimensionalities of contours. We see that, contrary to 
the expectations, the scalar isolines are smooth below L, 
which is quite natural from the physical perspective since 
all the factors (velocity, pumping, diffusion) are smooth 
at these scales. In particular, that means that the scalar 
field non-smoothness is related to the discontinuities of 
8(r) across (smooth) isolines. 

Let us now discuss the probability density function 
(PDF) of contour perimeters P and sizes characterized by 
the mean radius R = \J ((p — (p)) 2 ) , here p denotes the 
pairs of points parameterizing a contour and averaging (•) 
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FIG. 3. The generalized box counting fractal dimension D q 
is scale dependent. Below the forcing scale L the curves are 
smooth with D q — 1, while above L, the contour seems to be 
a simple fractal with D q in between 1.55 4- 1.7. The pumping 
frequency was v/X — 0.01. 
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FIG. 4. Box counting fractal dimension Do, for v/\ = 0.01, 
estimated from the derivative d log N(e)/ dlog(L/£). It is 
scale dependent and we see that bellow the forcing scale L the 
curves are smooth with Do — 1, while above L, the contour 
seems to be a fractal with Do in between 1.55 4- 1.7. 



is done over the points. Figures l5lrJ1 present the PDFs of 
\og(P/L) and \og(R/L) and show that they depend nei- 
ther on the resolution for rd/L = 0.015 -7- 0.06 nor on the 
pumping frequency for v j\ = 0.005 4- 0.05. In both fig- 
ures, the upper three curves are for different v /A and the 
same rd/L = 0.06; the lower three curves (PDF of those 
was divided by 10) differ in rd/L, while v/\ = 0.01. The 
only difference one can distinguish in the lower curves in 
Figures l5l6l is the appearance of the secondary maximum 
at small scales. This shows an abundance of diffusion- 
scale contours since the maximum appears only for the 
best resolution with the scale 0.015- L com parable to the 
molecular diffusion scale, which is -J Kd/ A = 0.0 IL in all 
runs. PDFs do not seem to depend on the velocity re- 
alization so Figures 15161 were obtained by averaging over 
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FIG. 5. The PDF of perimeters for different pumping frequen- 
cies (upper three curves) and resolutions (lower three curves 
shifted down by dividing by 10). 
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FIG. 6. The PDF of sizes for different pumping frequen- 
cies (upper three curves) and resolutions (lower three curves 
shifted down by dividing by 10). 



different realizations. 

Let us consider separately left and right tails of the 
PDFs. Since the probability of contours much larger than 
Td is independent of diffusion then it is determined by an 
interplay of stretching and pumping. Figures 15 161 show 
that the left tails of both V{P) and V{R) look like a 
power law with the power 3/2. Contours shorter than 
L must appear when pumping cuts a piece off a thin 
long contour, the probability of such a cut is oc P oc R 
(small contours are smooth so that P oc R). Extra factor 
VP oc VP in the PDF may appear because to be ob- 
served small contours need to survive without being swal- 
lowed by further pumping events, the lifetime is likely to 
be oc VP oc VP- Since creation and survival are inde- 
pendent events, their probabilities are multiplied. 

What one may expect for the PDFs of contour sizes 
at scales exceeding LI It is tempting to assume that a 
long contour appears due to an evolution undisturbed by 
pumping during a long time t, the length of such contour 
is Lexp(Ai) as long as it does not exceed L 2 /rd- The 
probability that pumping did not act during the time t 
is given by a Poisson law exp(— vt) where v is the pump- 
ing correlation time (in our algorithm, an average time 
between producing blobs of 9 in a given place). We then 
obtain V(R) = J dtcxp(~vt)S(R- exp(At)) oc R- X ~ u / X , 
which would mean that the PDF tail is non-universal 



and depends on the statistics of pumping and velocity. 
If that was true, the same tail one would expect for the 
PDF of perimeters as well. However, the above consid- 
eration totally disregards the fractal nature of long con- 
tours (shown in Figure [3]). We have checked that the 
fractal dimensions within our error bars were the same 
for different pumping frequencies v, sec Table [T] In line 
with long contour fractality, the right tails of the PDFs 
of P and R are very different. The tail of the PDF of 
sizes, V{R), looks log-normal, see Figure [5] The tail of 
V{P) looks like a universal power law, independent of the 
resolution and the pumping frequency v, see Figures [5] 
In log coordinates the tail is close to 1/P so it may be 
that V(P) oc P~ 2 , yet we were unable to derive this law 
theoretically so far. 
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0.05 


1.61 ±0.05 


1.66 ± 0.02 


1.67 ±0.02 


0.02 


1.62 ± 0.03 


1.65 ± 0.05 


1.65 ±0.05 


0.01 


1.62 ±0.05 


1.67 ±0.05 


1.67 ±0.03 


0.005 


1.62 ± 0.03 


1.68 ± 0.06 


1.68 ±0.04 



TABLE I. Generalized fractal dimensions, for different pump- 
ing frequencies u/\. 

Within our accuracy, we cannot see any difference be- 
tween the dimensionalities of the different orders and con- 
clude that our contours are mono-fractals in distinction 
from multi-fractal iso-vorticity contours in a direct cas- 
cade of 2d turbulence [H . This difference might be due to 
the fact that all parts of our scalar contours go through 
the same history of velocity, while parts of a long vortic- 
ity contour may have different histories. 

Looking at Figures 15161 a natural question is whether 
the positions of the maxima depend on the resolution. 
Pumping produces more or less circular contours of the 
radius L, which are then deformed into ellipsoids of in- 
creasing eccentricity by the velocity field. Pumping con- 
tinues to act by bending elongated contours. Those con- 
tours that on average have not changed much by this 
bending disappear after reaching the length of order 
L 2 /rd and respectively the width of order the coarse- 
graining scale rd- One then asks if the scale L 2 /rd is 
special, apart from L and r^. Numerics give a negative 
answer: lower curves in Figures 15161 compare runs with 
three different and show that the PDFs do not depend 
on the resolution for the scales exceeding r^. In partic- 
ular, the PDFs of the sizes have the maximum at the 
pumping scale L independently of the resolution (or dif- 
fusion scale) . We checked that the statistics is practically 
the same for either an average distance of the contour 
points from their center of mass or the maximal distance 
between the points of the contour (gyration radius). The 
PDFs of perimeter also have maxima independent of the 
diffusion scale (at the size approximately 2irL). The only 
difference one can distinguish in the lower curves in Fig- 
ures 15161 is the appearance of the secondary maximum at 
small scales for the curve with the best resolution, there 
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FIG. 7. Effective diffusivity of the driving function for 
velocity realizations (1-^3), where t is time in Lowner's eq., see 
[ij. The curves are ensemble averages of Ni contours (inset: 
of rii contours) from the i— th velocity realization. Inset shows 
the effective diffusivity for velocity realizations 1 and 3 after 
contraction by L /ra- 



the resolution scale 0.015L is comparable with the dif- 
fusion scale, which is 0.01L in all runs; in other words, 
there is an abundance of diffusion-scale contours. 

Quite unusual statistics of the passive scalar lacking 
scale-invariance has been found at r > L: multi-point 
correlation functions strongly depend on geometry [l8j]. 
This is because of the highly anisotropic " strip" structure 
of the field seen clearly in Figure [T] Let us discuss the 
statistics of bending for the isolines extending for such 
long distances. If there was only pumping then on the 
scales larger than L the scalar would be a short-correlated 
field whose nodal lines are equivalent to critical percola- 
tion i.e. SLEg. Without diffusion and with infinite res- 
olution, velocity only distorts the field. Of course, dis- 
torted field is not SLE 0]; f° r instance, stretching verti- 
cally a chordal SLE in a half-plane one adds to the driving 
function extra intervals of no change, that diminishes k 
and provides for a finite correlation scale (equivalcntly, 
finite correlations in Lowner time t which is the coordi- 
nate along the contour) . And yet deforming it back (with 
a time-dependent distortion factor cxp(2AT)) we would 
get the same SLE 6 . While the restoration procedure it- 
self may be not very practical since real flows consist of 
many such domains oriented randomly, the very possibil- 
ity of it means potential availability of very useful exact 
formulae describing the statistics of contours. For exam- 
ple, one may be able to describe crossing and surround- 
ing probabilities (like Cardy-Smirnov formula |20| ) 
with just a simple re-scaling. However, a finite resolu- 
tion/diffusion leads to irreversible changes (and makes 
the statistics stationary). Indeed, when velocity distorts 
the contours, it causes some distances in the contract- 
ing directions to become less than the resolution scale 
which leads to reconnections and disappearance of thin 
contours. Now it is highly nontrivial if any trace of initial 
SLE can be recovered, and if yes, after what contraction. 





N 1= 486, TX=100, 


0.8 




'— T /L = 0.06 

d 




My &fta*K 


N = 451, Tk = 100, 


0.6 




r /L = 0.03 

d 


A 




N 3 = 363, TX, = 100, 


^0.4 




r/L = 0.015 

d 






N 4 = 175, TX, = 50, 


0.21 




-e-r A = 0.06 

d 






N 5 = 88, TX = 20, 


» 


?Qm %m ~ - — - « . 


^r\/L = 0.06 

d 


c 


) 2000 4000 6000 8000 10000 





t[L 2 ] 



FIG. 8. Diffusivity k for different resolution scales and the 
evolution times TX for the same velocity realization. The val- 
ues Ni show the number of contours we averaged over and 
label the curves. Graphs Ni,N% and iVg show that k de- 
pends on the evolution time TX G {20,50,100}, while Ni, 
N2 and N3 show that k does not depend on the resolution 
Td/L £ {0.015, 0.03, 0.06}. The curves of different resolutions 
are of different length since the time of the Lowner map is 
measured in the units of length squared and for a better reso- 
lution we used a smaller physical window (all scalar snapshots 
had 40000 2 pixels, but in terms of the forcing scale they had 
(600£) 2 , (1200L) 2 , (2400L) 2 respectively). 



One may try to estimate the stretching factor empirically 
by measuring the aspect ratio of the boxes one can fit the 
contours into; we found out that this approach does not 
work since such an aspect ratio fluctuates strongly (by 
several orders of magnitude). 

Let us remind that Oseledec theorem and a general 
theory of Lagrangian chaos with diffusion state that 
any finite-point statistics is independent of the velocity 
realization and can be understood under the assumption 
that the effective contraction factor is L/r^. This is based 
on the fact that Lyapunov exponents are self-averaging 
and the mean lifetime of scalar blobs is A" 1 hx{L/rd). If 
one is to extend this reasoning to the statistics of an 
infinite-point object, then one expects the isoline statis- 
tics to depend strongly on the resolution scale and be 
independent of the velocity realization in a steady state 
and for sufficiently long contours that appear after a long 
history of stretching. As we now see, both predictions fail 
spectacularly. 

It was already stated above - we characterize the be- 
havior of a single line by the driving function of the 
Lowner map. Figure [7] shows that driving functions be- 
have roughly as that of diffusion and so can be char- 
acterized by the diffusivity k. It is independent of the 
resolution scale (this may be related to the problem of 
noise sensitivity in the statistics of random curves [2l|). 
as seen in Figure [51 It is also independent of the contour 
size as long as it is larger than L, yet k is different in 
different velocity realizations even on the same very long 
Lowner timescale (i.e. for very long contours). Most 
dramatically, Figure [5] shows that k depends strongly on 
time T (of stretching) on extremely long timescalcs far 
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exceeding the time A" 1 \n(L/rd) during which the scalar 
field itself acquires stationary statistics. We conclude 
that long contours undergo stretching and their statis- 
tics changes on a much larger timescale than the scalar 
blob lifetime. Let us stress the discrepancy: k depends 
strongly on the distortion and yet it is independent of rn 
(at least within the limits we studied), that is the effec- 
tive distortion factor for long isolines is not L/rd (as it is 
for multi-point statistics). On the other hand, note that 
choosing the contraction factor equal to L/r d one ob- 
tains k of the "right" order of magnitude (between 3 and 
17 for different velocity realizations, see inset of Figure 
[7]). One is tempted to explore whether one can recover 



SLE contours (despite all the loss of information due to 
reconnections) by fine-tuning the distortion factor, possi- 
bly by requiring either restoration of statistical isotropy 
or shortest correlation time of £(i). Further studies with 
extensive statistics are needed to sort out which proper- 
ties of the critical percolation are retained by the large- 
scale statistics of passive scalar contours in the Batchelor 
regime. 
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